{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\"> \n",
    "</center>\n",
    "\n",
    "## <center> Открытый курс по машинному обучению </center>\n",
    "\n",
    "### <center> Автор материала: Тетерников Илья (@tetelias) </center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center>Индивидуальный проект по анализу данных:</center>\n",
    "## <center>Лихорадка Денгэ: предсказание распространения</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Описание набора данных и признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Введение\n",
    "\n",
    "Данный материал описывает участие в соревновании \"DengAI: Predicting Disease Spread\" на портале Drivedata.com. \n",
    "Для исследования предоставляется набор данных о выявленных случаях заражения в двух городах: Сан-Хуан, Пуэрто-Рико и Икитос, Перу в течение 18 и 10 лет соответственно. Население городов составляет 395 и 406 тыс. человек, т.е. количества случаев заболевания можно не нормализовывать.\n",
    "\n",
    "Целью исследования является исследование зависимостей всплесков частоты заболевания от погодных условий. Для проверки требуется предсказать количество выявленных заболеваний в двух городах за период в 5 и 3 года соответственно.\n",
    "\n",
    "Информация о заболевании из [Википедии](https://en.wikipedia.org/wiki/Dengue_fever):\n",
    "\n",
    "* Частота заболевания увеличилась в 30 раз за последние 50 лет, и ежегодный ущерб от болезни исчисляется тысячами жизней и миллионами болезней\n",
    "* Территория распространения болезни расширяется за счет умеренных широт\n",
    "* Инкубационный период составляет в среднем 1 неделю\n",
    "* Болезнь легко распространяется и потому может приводить к эпидемиям\n",
    "\n",
    "Данный вирус переносится комарами рода Aedes, чаще всего вида Aedes aegypti. По этой причине динамика переноса лихорадки связана с климатическими условиями, такими как температура и количество осадков. \n",
    "\n",
    "Информация о переносчике заболевания из [Википедии](https://en.wikipedia.org/wiki/Aedes_aegypti): \n",
    "\n",
    "* является переносчиком этого и других опасных инфекционных заболеваний\n",
    "* жизненный цикл - 2-4 недели\n",
    "* проявляет активность в основном на рассвете и закате \n",
    "* проявляет активность круглый год \n",
    "* предпочитает стоячую воду\n",
    "* увеличение численности насекомых связывают с увеличением урбанизации\n",
    "\n",
    "\n",
    "Данные представляют собой объединение информации из различных источников: Центр по Контролю и Предотврщению Заболеваний США(CDC USA) и Национальное управление океанических и атмосферных исследований США(NOAA).\n",
    "\n",
    "Глубокое понимание зависимости болезни от погодных условий должно помочь в предотвращении эпидемий. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Целевой признак\n",
    "\n",
    "__total_cases__ - количество выявленных за соответствующую неделю заболеваний\n",
    "    \n",
    "__Тестовая выборка__ - последовательная по времени выборка из будущего по отношению к доступным для обучения данным."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Признаки:\n",
    "\n",
    "__city__ – название города: sj - Сан-Хуан и iq - Икитос  \n",
    "__year__ – год  \n",
    "__weekofyear__ – номер недели в текущем году  \n",
    "__week_start_date__ – дата в формате гггг-мм-дд\n",
    "\n",
    "##### Дневные замеры погодных условий в понедельник каждой недели:\n",
    "\n",
    "__station_max_temp_c__ – максимальная температура воздуха в течение суток в градусах Цельсия  \n",
    "__station_min_temp_c__ – минимальная температура воздуха в течение суток в градусах Цельсия  \n",
    "__station_avg_temp_c__ – средняя температура воздуха в течение суток в градусах Цельсия  \n",
    "__station_precip_mm__ – размер осадков, измеренный датчиками на поверхности  \n",
    "__station_diur_temp_rng_c__ – дневные колебания температуры, т.е. изменение температуры в период от восхода до заката, в градусах Цельсия  \n",
    "\n",
    "__precipitation_amt_mm__ – замер осадков, но произведенный с помощью спутника  \n",
    "\n",
    "##### Прогнозы погоды:\n",
    "\n",
    "__reanalysis_sat_precip_amt_mm__ – количество осадков в мм  \n",
    "__reanalysis_dew_point_temp_k__ – средняя температура росы в градусах Кельвина  \n",
    "__reanalysis_air_temp_k__ – средняя температура воздуха в течение суток на уровне 2м от земли в градусах Кельвина  \n",
    "__reanalysis_relative_humidity_percent__ – относительная влажность  \n",
    "__reanalysis_specific_humidity_g_per_kg__ – абсолютная влажность в г воды на кг воздуха  \n",
    "__reanalysis_precip_amt_kg_per_m2__ – количество осадков в кг на м2  \n",
    "__reanalysis_max_air_temp_k__ – максимальная температура воздуха в течение суток в градусах Кельвина  \n",
    "__reanalysis_min_air_temp_k__ – минимальная температура воздуха в течение суток в градусах Кельвина  \n",
    "__reanalysis_avg_temp_k__ – средняя температура воздуха в течение суток на уровне земл в градусах Кельвина  \n",
    "__reanalysis_tdtr_k__ – дневные колебания температуры в градусах Кельвина  \n",
    "\n",
    "##### Нормализованный ровень растительности на территории:\n",
    "\n",
    "__ndvi_se__ – юго-восточной части города  \n",
    "__ndvi_sw__ – юго-западной части города  \n",
    "__ndvi_ne__ – северо-восточной части города  \n",
    "__ndvi_nw__ – северо-западной части города  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Первичный анализ данных "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Загрузка данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Импортируем необходимые библиотеки для обработки и визуализации данных\n",
    "from hyperopt import fmin, tpe, rand, hp, STATUS_OK, Trials\n",
    "from mlxtend.feature_selection import SequentialFeatureSelector as SFS\n",
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.model_selection import TimeSeriesSplit, cross_val_score, GridSearchCV, learning_curve\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from warnings import filterwarnings\n",
    "filterwarnings('ignore')\n",
    "\n",
    "import math\n",
    "import matplotlib.pylab as pylab\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import xgboost as xgb\n",
    "plt.style.use('ggplot')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные предоставлены в виде 3 файлов: признаки для обучения, целевой признак и признаки для предсказания."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = pd.read_csv('../../data/dengue_features_train.csv', infer_datetime_format=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = pd.read_csv('../../data/dengue_labels_train.csv', infer_datetime_format=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = pd.read_csv('../../data/dengue_features_test.csv', infer_datetime_format=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df = pd.merge(features, labels, how='outer', on=labels.columns.tolist()[:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df = test.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df.iloc[:, :15].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df.iloc[:, 15:].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Исследование признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Два города сильно различаются по климату и качеству сбора информации, поэтому будем рассматривать признаки отдельно по городам"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверим, все ли данные в наличии:\n",
    "\n",
    "* как видим, в Сан-Хуане много пропущенных значений по растительности, а в Икитос отсутствует значительная часть наземных замеров климатических условий. \n",
    "* есть проблемы и с другими замерами, но там это составляет меньше 1% от общего числа.\n",
    "* в тестовых выборках ситуация гораздо лучше: проблемы по городам те же самые, но гораздо менее выраженные\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='sj'].isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='iq'].isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df[test_df.city=='sj'].isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df[test_df.city=='iq'].isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Взглянем на переходы из года в год:\n",
    "* становится понятна природа части недостающих данных: люди празднуют новый год и им не до фиксирования замеров\n",
    "* некоторые недели приписаны к следующему году, что нужно исправить\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[(train_df.weekofyear==1)|(train_df.weekofyear==52)|(train_df.weekofyear==53)].iloc[:, :15]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df.year = np.where((train_df.week_start_date.apply(lambda x: x[5:]) == '01-01')&(train_df.weekofyear != 1),\n",
    "                         train_df.year-1, train_df.year)\n",
    "test_df.year = np.where((test_df.week_start_date.apply(lambda x: x[5:]) == '01-01')&(test_df.weekofyear != 1),\n",
    "                         test_df.year-1, test_df.year)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df = train_df.sort_values(by='week_start_date')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df.week_start_date = train_df.week_start_date.apply(pd.to_datetime)\n",
    "train_df.set_index('week_start_date', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df.set_index('week_start_date', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь посмотрим на распределение признаков:\n",
    "\n",
    "* значения признаков не противоречат природе данных\n",
    "* уровень растительности в Сан-Хуане принимает отрицательные значения, что связано с нормализацией\n",
    "* разброс температур и количество осадков значительно меньше в Сан-Хуане\n",
    "* в признаках, связанных с растительностью в Сан-Хуане и количеством осадков в обоих городах, мы наблюдаем выбросы\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = train_df.columns[1:3].tolist() + train_df.columns[3:25].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='sj'][cols[:13]].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='iq'][cols[:13]].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='sj'][cols[13:]].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='iq'][cols[13:]].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Найдем пары признаков с максимальной корреляцией:\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mostcorr(dataframe, numtoreport):\n",
    "    \n",
    "    cormatrix = dataframe.corr()\n",
    "    \n",
    "    cormatrix *= np.tri(*cormatrix.values.shape, k=-1).T\n",
    "    \n",
    "    cormatrix = cormatrix.stack()\n",
    "    cormatrix = cormatrix.reindex(cormatrix.abs().sort_values(ascending=False).index).reset_index()\n",
    "    \n",
    "    cormatrix.columns = [\"Признак №1\", \"Признак №2\", \"Корреляция\"]\n",
    "    return cormatrix.head(numtoreport)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mostcorr(train_df[train_df.city=='sj'],10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mostcorr(train_df[train_df.city=='iq'],10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Температура в Сан-Хуане изменяется слабо в течение дня, а абсолютная влажность движется вслед за температурой. \n",
    "* В Икитос количество вегетации в различных частях города имеет очень сильную корреляцию.\n",
    "* Наконец, нельзя не отметить, что три пары с наибольше корреляцией совпадают для обоих городов и едва отличимы от 1, так что удалим reanalysis_sat_precip_amt_mm, reanalysis_dew_point_temp_k, reanalysis_avg_temp_k:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del train_df['reanalysis_sat_precip_amt_mm']\n",
    "del train_df['reanalysis_dew_point_temp_k']\n",
    "del train_df['reanalysis_avg_temp_k']\n",
    "del test_df['reanalysis_sat_precip_amt_mm']\n",
    "del test_df['reanalysis_dew_point_temp_k']\n",
    "del test_df['reanalysis_avg_temp_k']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Посмотрим на связь признаков с целевой переменной:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale = 1.5)\n",
    "(abs(train_df[(train_df.city=='sj')].corr())\n",
    " .total_cases\n",
    " .drop('total_cases')\n",
    " .sort_values()\n",
    " .plot\n",
    " .barh())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale = 1.5)\n",
    "(abs(train_df[(train_df.city=='iq')].corr())\n",
    " .total_cases\n",
    " .drop('total_cases')\n",
    " .sort_values()\n",
    " .plot\n",
    " .barh())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Часть признаков присутствует в обоих группах.\n",
    "Это означает, что абсолютная влажность и оба замера минимальной температуры воздуха могут быть наиболее полезны в предсказании целевого признака. Интересно также присутствие признака 'year'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Исследование целевого признака"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим график относительной частотности количества случаев заболевания для двух городов. Оба распределения явно относятся к отрицательным биномиальным."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.FacetGrid(train_df, hue='city', aspect=4, size=4) \n",
    "fig.map(sns.kdeplot,'total_cases',shade=True)\n",
    "max_x = train_df.total_cases.max()\n",
    "min_x = train_df.total_cases.min()\n",
    "fig.set(xlim=(min_x,max_x))\n",
    "fig.set(ylim=(0, 0.09))\n",
    "fig.add_legend()\n",
    "fig.fig.suptitle(\"Количество случаев\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Несмотря на визуальную разницу распределения отличаются только масштабом:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df['total_cases_temp'] = train_df.total_cases.where(train_df.city=='sj', train_df.total_cases*3.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.FacetGrid(train_df, hue='city', aspect=4, size=4) \n",
    "fig.map(sns.kdeplot,'total_cases_temp',shade=True)\n",
    "max_x = train_df.total_cases.max()\n",
    "min_x = train_df.total_cases.min()\n",
    "fig.set(xlim=(min_x,max_x))\n",
    "fig.add_legend()\n",
    "fig.fig.suptitle(\"Количество случаев масштабированное\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del train_df['total_cases_temp']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вычислим параметры распределений для целевой переменной: скошенность, а также p и r."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_std = train_df[train_df.city=='sj'].total_cases.std()\n",
    "sj_mean = train_df[train_df.city=='sj'].total_cases.mean()\n",
    "sj_median = train_df[train_df.city=='sj'].total_cases.median()\n",
    "iq_std = train_df[train_df.city=='iq'].total_cases.std()\n",
    "iq_mean = train_df[train_df.city=='iq'].total_cases.mean()\n",
    "iq_median = train_df[train_df.city=='iq'].total_cases.median()\n",
    "\n",
    "#Скошенность\n",
    "sj_skew = (sj_mean - sj_median) / sj_std\n",
    "iq_skew = (iq_mean - iq_median) / iq_std\n",
    "\n",
    "#p\n",
    "sj_p = 1 - sj_mean / sj_std**2\n",
    "iq_p = 1 - iq_mean / iq_std**2\n",
    "\n",
    "#r\n",
    "sj_r = sj_mean**2 / sj_std**2 / sj_p\n",
    "iq_r = iq_mean**2 / iq_std**2 / iq_p\n",
    "\n",
    "print('Скошенность распределения целевой переменной для Сан-Хуана равна:%0.3f' % sj_skew)\n",
    "print('Скошенность распределения целевой переменной для Икитос равна:%0.3f' % iq_skew)\n",
    "print('')\n",
    "print('Параметры p и r распределения целевой переменной для Сан-Хуана равны:%0.3f' % sj_p + ' и %0.3f' % sj_r)\n",
    "print('Параметры p и r распределения целевой переменной для Икитос равны:%0.3f' % iq_p + ' и %0.3f' % iq_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Это еще одно доказательство схожести распределений целевой переменной для обоих городов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на историческое распределение случаев заболевания в обоих городах:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='sj'].plot(kind='line',  x=['year', 'weekofyear'], y='total_cases',\n",
    "              figsize=[15,7],title='Распределение случаев болезни в Сан-Хуане по времени')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df[train_df.city=='iq'].plot(kind='line',  x=['year', 'weekofyear'], y='total_cases',\n",
    "              figsize=[15,7],title='Распределение случаев болезни в Икитос по времени')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные по Икитос шумнее, но в целом видно, что всплески заболевания имеют схожую продолжительность."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь посмотрим на распределение случаев заболевания понедельно:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.FacetGrid(train_df, hue='city', aspect=4, size=4) \n",
    "fig.map(sns.pointplot,'weekofyear','total_cases')\n",
    "max_x = train_df.weekofyear.max()\n",
    "min_x = train_df.weekofyear.min()\n",
    "fig.set(xlim=(min_x,max_x))\n",
    "fig.set(ylim=(0, 130))\n",
    "fig.add_legend()\n",
    "fig.fig.suptitle(\"Распределение случаев болезни в течение года\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наблюдается интересная сезонность с минимумом в районе 13-19 недель и максимумом - 40-46. Более точные выводы невозможно сделать ввиду малого размера выборки."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для приближения распределения целевого признака к нормальной форме используем операцию логарифмирования:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df['total_cases_adj'] = train_df.total_cases.apply(math.log1p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Первичный визуальный анализ данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Принимая во внимание срок инкубации и жизненный цикл распространителя, попробуем рассмотреть влияние текущих погодных условий на уровень заболевания в будущем на 1, 2, 3 и 4 недели. Для этого создадим новые переменные из целевой сдвигом на соответствующее число недель."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_df = train_df[train_df.city=='sj'].copy()\n",
    "iq_df = train_df[train_df.city=='iq'].copy()\n",
    "for i in range(1, 5):\n",
    "    sj_df['total_cases'+str(i)] = sj_df.total_cases.shift(i)\n",
    "    iq_df['total_cases'+str(i)] = iq_df.total_cases.shift(i)\n",
    "    \n",
    "train_df = pd.concat((sj_df, iq_df))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так как города имеют разный климат, взглянем на взаимодействие признаков для каждого из них отдельно:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(30, 30))\n",
    "corr = train_df[(train_df.city=='sj')].corr()\n",
    "sns.heatmap(corr, mask=np.zeros_like(corr, dtype=np.bool), \n",
    "            cmap=sns.diverging_palette(220, 10, as_cmap=True), annot=True, square=True, center=0, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(30, 30))\n",
    "corr = train_df[train_df.city=='iq'].corr()\n",
    "sns.heatmap(corr, mask=np.zeros_like(corr, dtype=np.bool), \n",
    "            cmap=sns.diverging_palette(220, 10, as_cmap=True), annot=True, center=0, square=True, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Удалим созданные сдвиговые переменные:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del train_df['total_cases1']\n",
    "del train_df['total_cases2']\n",
    "del train_df['total_cases3']\n",
    "del train_df['total_cases4']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Посмотрим теперь на распределение признаков:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc(\"font\", size=13)\n",
    "plt.figure(figsize=(20,50))\n",
    "alpha=0.6\n",
    "\n",
    "\n",
    "ax1 = plt.subplot2grid((6,3), (0,0))\n",
    "train_df.ndvi_ne.plot(kind=\"kde\",alpha=alpha)\n",
    "ax1.set_xlabel(\"ndvi_ne\")\n",
    "\n",
    "ax2 = plt.subplot2grid((6,3),(0,1))\n",
    "train_df.ndvi_nw.plot(kind=\"kde\",alpha=alpha)\n",
    "ax2.set_xlabel(\"ndvi_nw\")\n",
    "\n",
    "ax3 = plt.subplot2grid((6,3),(0,2))\n",
    "train_df.ndvi_se.plot(kind=\"kde\",alpha=alpha)\n",
    "ax3.set_xlabel(\"ndvi_se\")\n",
    "\n",
    "ax4 = plt.subplot2grid((6,3),(1,0))\n",
    "train_df.ndvi_sw.plot(kind=\"kde\",alpha=alpha)\n",
    "ax4.set_xlabel(\"ndvi_sw\")\n",
    "\n",
    "ax5 = plt.subplot2grid((6,3),(1,1))\n",
    "train_df.reanalysis_air_temp_k.plot(kind=\"kde\",alpha=alpha)\n",
    "ax5.set_xlabel(\"reanalysis_air_temp_k\")\n",
    "\n",
    "ax6 = plt.subplot2grid((6,3),(1,2))\n",
    "train_df.reanalysis_max_air_temp_k.plot(kind=\"kde\",alpha=alpha)\n",
    "ax6.set_xlabel(\"reanalysis_max_air_temp_k\")\n",
    "\n",
    "ax7 = plt.subplot2grid((6,3),(2,0))\n",
    "train_df.reanalysis_min_air_temp_k.plot(kind=\"kde\",alpha=alpha)\n",
    "ax7.set_xlabel(\"reanalysis_min_air_temp_k\")\n",
    "\n",
    "ax8 = plt.subplot2grid((6,3),(2,1))\n",
    "train_df.reanalysis_precip_amt_kg_per_m2.plot(kind=\"kde\",alpha=alpha)\n",
    "ax8.set_xlabel(\"reanalysis_precip_amt_kg_per_m2\")\n",
    "\n",
    "ax9 = plt.subplot2grid((6,3),(2,2))\n",
    "train_df.reanalysis_relative_humidity_percent.plot(kind=\"kde\",alpha=alpha)\n",
    "ax9.set_xlabel(\"reanalysis_relative_humidity_percent\")\n",
    "\n",
    "ax10 = plt.subplot2grid((6,3),(3,0))\n",
    "train_df.precipitation_amt_mm.plot(kind=\"kde\",alpha=alpha)\n",
    "ax10.set_xlabel(\"precipitation_amt_mm\")\n",
    "\n",
    "ax11 = plt.subplot2grid((6,3),(3,1))\n",
    "train_df.reanalysis_specific_humidity_g_per_kg.plot(kind=\"kde\",alpha=alpha)\n",
    "ax11.set_xlabel(\"reanalysis_specific_humidity_g_per_kg\")\n",
    "\n",
    "ax12 = plt.subplot2grid((6,3),(3,2))\n",
    "train_df.reanalysis_tdtr_k.plot(kind=\"kde\",alpha=alpha)\n",
    "ax12.set_xlabel(\"reanalysis_tdtr_k\")\n",
    "\n",
    "ax13 = plt.subplot2grid((6,3),(4,0))\n",
    "train_df.station_avg_temp_c.plot(kind=\"kde\",alpha=alpha)\n",
    "ax13.set_xlabel(\"station_avg_temp_c\")\n",
    "\n",
    "ax14 = plt.subplot2grid((6,3),(4,1))\n",
    "train_df.station_diur_temp_rng_c.plot(kind=\"kde\",alpha=alpha)\n",
    "ax14.set_xlabel(\"station_diur_temp_rng_c\")\n",
    "\n",
    "ax15 = plt.subplot2grid((6,3),(4,2))\n",
    "train_df.station_max_temp_c.plot(kind=\"kde\",alpha=alpha)\n",
    "ax15.set_xlabel(\"station_max_temp_c\")\n",
    "\n",
    "ax16 = plt.subplot2grid((6,3),(5,0))\n",
    "train_df.station_min_temp_c.plot(kind=\"kde\",alpha=alpha)\n",
    "ax16.set_xlabel(\"station_min_temp_c\")\n",
    "\n",
    "ax17 = plt.subplot2grid((6,3),(5,1))\n",
    "train_df.station_precip_mm.plot(kind=\"kde\",alpha=alpha)\n",
    "ax17.set_xlabel(\"station_precip_mm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "По своей природе признаки делятся на три части: уровень вегетации, температура, влажность и количество осадков. Для первых трех наблюдается распределение близкое к нормальному, а вот последняя группа имеет биномиальное распределение.\n",
    "\n",
    "Воспользуемся логарифмированием для нормализации их распределения:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_df.precipitation_amt_mm = train_df.precipitation_amt_mm.apply(math.log1p)\n",
    "test_df.precipitation_amt_mm = test_df.precipitation_amt_mm.apply(math.log1p)\n",
    "train_df.station_precip_mm = train_df.station_precip_mm.apply(math.log1p)\n",
    "test_df.station_precip_mm = test_df.station_precip_mm.apply(math.log1p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Инсайты, найденные зависимости"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Ни один из признаков не имеет сильной связи с целевой переменной\n",
    "* Признаки температуры, что следует из природы признаков, сильно коррелируют"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* учитывая разницу в расположении относительно моря, разброс температур и количество осадков значительно меньше в Сан-Хуане\n",
    "* Наличие морского климата приводят к тому, что температура в Сан-Хуане изменяется слабо в течение дня, а абсолютная влажность движется вслед за температурой.\n",
    "* а вот максимальное количество осадков больше все же в столице Пуэрто-Рико из-за расположения в зоне торнадо. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Визуальный анализ:\n",
    "\n",
    "* нормализованное распределение целевой переменной гораздо лучше коррелирует с большинством признаков, чем скошенное\n",
    "* для Икитос четко прослеживается повышение корреляции признаков с целевой переменной при смещении на 1, 2 и 3 недели и небольшое падение на последнем сдвиге. А вот для Сан-Хуана отсутствие сдвига дает самую лучший набор корреляций целевой переменной с признаками. Так что для Икитос при построении модели имеет смысл сдвинуть наборы переменных на 3 недели назад.\n",
    "* Сложно сказать, с чем связана закономерность из предыдущего пункта. Возможно, \n",
    "* не видно заметной зависимости числа заболеваний от уровня растительности. Учитывая информацию из Википедии, что основной источник размножения переносчиков заболевания - это стоячие водоемы вблизи/в жилых районах, данные признаки представляются несущественными и их можно удалить.\n",
    "* Из двух признаков, характеризующих дневные колебания, reanalysis_tdtr_k показывает лучшую корреляцию. Удалим station_diur_temp_rng_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del train_df['ndvi_ne']\n",
    "del train_df['ndvi_nw']\n",
    "del train_df['ndvi_sw']\n",
    "del train_df['ndvi_se']\n",
    "del train_df['station_diur_temp_rng_c']\n",
    "\n",
    "del test_df['ndvi_ne']\n",
    "del test_df['ndvi_nw']\n",
    "del test_df['ndvi_sw']\n",
    "del test_df['ndvi_se']\n",
    "del test_df['station_diur_temp_rng_c']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Выбор метрики"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выбор метрики задан организаторами: средняя абсолютная ошибка(mean absolute error, MAE).\n",
    "Она вычисляется по формуле:\n",
    "\n",
    "$$\\large \\begin{array}{rcl}\\mathcal{MAE} &=& \\frac{1}{n} \\sum_{i=1}^n \\left |f_i - y_i \\right|\n",
    "\\end{array}$$\n",
    "\n",
    "Обоснование выбора именно этой метрики организаторами дано не было, однако можно предположить, что они сделали выбор, основываясь на главном преимуществе MAE - простоте интерпретации.  \n",
    "\n",
    "Учитывая заявленное желание предсказывать эпидемии, использование средней квадратичной ошибки, которая гораздо лучше отслеживает выбросы, которыми являются вспышки эпидемий, было бы предпочтительнее."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Выбор модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В рамках данного проекта нам необходимо решить задачу регрессии.  \n",
    "Выбор ограничен сложными начальными условиями: небольшие набор данных и количество признаков.  \n",
    "\n",
    "Как потенциальные варианты решения рассмотрим:  \n",
    "* __Lasso Regression__ - регрессор c l1-регуляризацией  \n",
    "* __XGBRegressor__ - регрессор, использующий деревья, на основе алгоритма xgboost.  \n",
    "\n",
    "Посмотрим, какие результаты дадут модели на кросс-валидации и выберем лучшую."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Предобработка данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = ['year', 'weekofyear', 'precipitation_amt_mm', 'reanalysis_air_temp_k', 'reanalysis_max_air_temp_k',\n",
    "       'reanalysis_min_air_temp_k', 'reanalysis_precip_amt_kg_per_m2', 'reanalysis_relative_humidity_percent',\n",
    "       'reanalysis_specific_humidity_g_per_kg', 'reanalysis_tdtr_k', 'station_avg_temp_c', 'station_max_temp_c', \n",
    "        'station_min_temp_c', 'station_precip_mm', 'total_cases', 'total_cases_adj']\n",
    "\n",
    "features_to_scale = ['year', 'weekofyear', 'precipitation_amt_mm', 'reanalysis_air_temp_k', 'reanalysis_max_air_temp_k', \n",
    "                     'reanalysis_min_air_temp_k', 'reanalysis_precip_amt_kg_per_m2', \n",
    "                     'reanalysis_relative_humidity_percent', 'reanalysis_specific_humidity_g_per_kg', \n",
    "                     'reanalysis_tdtr_k', 'station_avg_temp_c', 'station_max_temp_c', \n",
    "                     'station_min_temp_c', 'station_precip_mm']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Разделим данные по городам:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_df = pd.concat([train_df[train_df.city=='sj'].copy(), test_df[test_df.city=='sj'].copy()])\n",
    "\n",
    "iq_df = pd.concat([train_df[train_df.city=='iq'].copy(), test_df[test_df.city=='iq'].copy()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_df = sj_df.drop(['city'], axis=1)\n",
    "iq_df = iq_df.drop(['city'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_sj = train_df[train_df.city=='sj'].shape[0]\n",
    "id_iq = train_df[train_df.city=='iq'].shape[0]\n",
    "vid_sj = int(0.7*id_sj)\n",
    "vid_iq = int(0.7*id_iq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сдвинем данные для Икитос на 3 недели назад. Образовавшийся пробел ввиду природы данных лучше заполнить последними значениями. Оставим также изначальный вариант, чтобы проверить гипотезу о сдвиге."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iq_base_df = iq_df.copy()\n",
    "iq_df[cols] = iq_df[cols].shift(-3)\n",
    "\n",
    "iq_df.iloc[iq_df.shape[0]-3] = iq_df.iloc[iq_df.shape[0]-4].values\n",
    "iq_df.iloc[iq_df.shape[0]-2] = iq_df.iloc[iq_df.shape[0]-4].values\n",
    "iq_df.iloc[iq_df.shape[0]-1] = iq_df.iloc[iq_df.shape[0]-4].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Заполним пропущенные данные методом заполнения вперед: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_df.fillna(method='ffill', inplace=True)\n",
    "iq_df.fillna(method='ffill', inplace=True)\n",
    "iq_base_df.fillna(method='ffill', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Отмасштабируем признаки:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "features_scaled_sj = sj_df[features_to_scale].values\n",
    "sjx_train = scaler.fit_transform(features_scaled_sj[:id_sj])\n",
    "sjx_test = scaler.transform(features_scaled_sj[id_sj:])\n",
    "sj_df[features_to_scale] = np.concatenate([sjx_train, sjx_test])\n",
    "\n",
    "features_scaled_iq = iq_df[features_to_scale].values\n",
    "iqx_train = scaler.fit_transform(features_scaled_iq[:id_iq])\n",
    "iqx_test = scaler.transform(features_scaled_iq[id_iq:])\n",
    "iq_df[features_to_scale] = np.concatenate([iqx_train, iqx_test])\n",
    "\n",
    "features_scaled_iq_base = iq_base_df[features_to_scale].values\n",
    "iq_base_train = scaler.fit_transform(features_scaled_iq_base[:id_iq])\n",
    "iq_base_test = scaler.transform(features_scaled_iq_base[id_iq:])\n",
    "iq_base_df[features_to_scale] = np.concatenate([iq_base_train, iq_base_test])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Учитывая, что мы имеем дело с временным рядом, корректным является только выделение валидационной выборки из конца обучающей:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sjx_valid_df = pd.DataFrame(sjx_train[vid_sj:id_sj], columns=features_to_scale)\n",
    "sjx_valid_df['y'] = sj_df.total_cases[vid_sj:id_sj].values\n",
    "sjx_valid_df['y_adj'] = sj_df.total_cases_adj[vid_sj:id_sj].values\n",
    "\n",
    "sjx_train_df = pd.DataFrame(sjx_train[:vid_sj], columns=features_to_scale)\n",
    "sjx_train_df['y'] = sj_df.total_cases[:vid_sj].values\n",
    "sjx_train_df['y_adj'] = sj_df.total_cases_adj[:vid_sj].values\n",
    "\n",
    "sjx_test_df = pd.DataFrame(sjx_test, columns=features_to_scale)\n",
    "\n",
    "iqx_valid_df = pd.DataFrame(iqx_train[vid_iq:id_iq], columns=features_to_scale)\n",
    "iqx_valid_df['y'] = iq_df.total_cases[vid_iq:id_iq].values\n",
    "iqx_valid_df['y_adj'] = iq_df.total_cases_adj[vid_iq:id_iq].values\n",
    "\n",
    "iqx_train_df = pd.DataFrame(iqx_train[:vid_iq], columns=features_to_scale)\n",
    "iqx_train_df['y'] = iq_df.total_cases[:vid_iq].values\n",
    "iqx_train_df['y_adj'] = iq_df.total_cases_adj[:vid_iq].values\n",
    "\n",
    "iqx_test_df = pd.DataFrame(iqx_test, columns=features_to_scale)\n",
    "\n",
    "iqx_base_valid_df = pd.DataFrame(iq_base_train[vid_iq:id_iq], columns=features_to_scale)\n",
    "iqx_base_valid_df['y'] = iq_base_df.total_cases[vid_iq:id_iq].values\n",
    "iqx_base_valid_df['y_adj'] = iq_base_df.total_cases_adj[vid_iq:id_iq].values\n",
    "\n",
    "iqx_base_train_df = pd.DataFrame(iq_base_train[:vid_iq], columns=features_to_scale)\n",
    "iqx_base_train_df['y'] = iq_base_df.total_cases[:vid_iq].values\n",
    "iqx_base_train_df['y_adj'] = iq_base_df.total_cases_adj[:vid_iq].values\n",
    "\n",
    "iqx_base_test_df = pd.DataFrame(iq_base_test, columns=features_to_scale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Кросс-валидация и настройка гиперпараметров модели (4 балла)\n",
    "\n",
    "    (+) Кросс-валидация выполнена технически верно, нет утечек данных. Разумно выбрано количество фолдов и разбиение (Random/Stratified или иное), зафиксирован seed. Присутствует объяснение. Объяснены гиперпараметры модели и способ их выбора. Выбор основан на некотором исследовании гипрепараметров модели для данной задачи;\n",
    "    \n",
    "    (+/-) Присутствуют незначительные ошибки (например, не зафиксирован seed) или отсутствует объяснение кросс-валидации. Но гиперпараметры модели и способ их выбора должны быть объяснены;\n",
    "    \n",
    "    (-/+) Кросс-валидация выполнена со значительными ошибками (например, преобразования данных проводится на всей выборке, таким образом возникает утечка данных из тестовой части выборки и, соответственно, результат кросс-валидации может иметь слишком оптимистичное значение). Гиперпараметры модели и способ их выбора не объяснены;\n",
    "    \n",
    "    (-) Отсутствуют."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выделим обучающую и тестовую выборки, а также целевую переменную для трех наборов данных:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_sj = sj_df.iloc[:id_sj].drop(['total_cases', 'total_cases_adj'], axis=1)\n",
    "y_sj = sj_df.iloc[:id_sj].total_cases_adj\n",
    "y_true_sj = sj_df['total_cases'].iloc[:id_sj]\n",
    "test_sj = sj_df.iloc[id_sj:].drop(['total_cases', 'total_cases_adj'], axis=1)\n",
    "\n",
    "X_iq = iq_df.iloc[:id_iq].drop(['total_cases', 'total_cases_adj'], axis=1)\n",
    "y_iq = iq_df.iloc[:id_iq].total_cases_adj\n",
    "y_true_iq = iq_df['total_cases'].iloc[:id_iq]\n",
    "test_iq = iq_df.iloc[id_iq:].drop(['total_cases', 'total_cases_adj'], axis=1)\n",
    "\n",
    "X_iq_base = iq_base_df.iloc[:id_iq].drop(['total_cases', 'total_cases_adj'], axis=1)\n",
    "y_iq_base = iq_base_df.iloc[:id_iq].total_cases_adj\n",
    "y_true_iq_base = iq_base_df['total_cases'].iloc[:id_iq]\n",
    "test_iq_base = iq_base_df.iloc[id_iq:].drop(['total_cases', 'total_cases_adj'], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Будем производить поиск гиперпараметров XGBRegressor с помощью GridSearchCV, используя методику TimeSeriesSplit для осуществления кросс-валидации."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Также будем осуществлять поиск лучшего набора признаков при помощи метода SFS библиотеки mlxtend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tss = TimeSeriesSplit(n_splits=10)\n",
    "set_params = {'n_estimators': [100, 200, 400],\n",
    "              'learning_rate': [0.1, 0.3, 0.5, 1],\n",
    "              'max_depth': [4, 5, 6]\n",
    "             }\n",
    "clf_xgb = xgb.XGBRegressor(objective='reg:linear', seed=13)\n",
    "grid = GridSearchCV(clf_xgb, \n",
    "                    set_params,\n",
    "                    cv=tss, \n",
    "                    scoring='neg_mean_absolute_error',\n",
    "                    n_jobs=-1\n",
    "                   )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.fit(X_sj, y_sj)\n",
    "sj_best = grid.best_estimator_\n",
    "pred_sj = sj_best.predict(test_sj)\n",
    "\n",
    "SFS1 = SFS(sj_best, k_features=(1, 14), forward=True, floating=True,\n",
    "           verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)\n",
    "sfs1 = SFS1.fit(X_sj.as_matrix(), y_sj)\n",
    "print([X_sg.columns[i-1] for i in sfs1.k_features])\n",
    "print(sfs1.k_score_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.fit(X_iq, y_iq)\n",
    "iq_best = grid.best_estimator_\n",
    "pred_iq = iq_best.predict(test_iq)\n",
    "print(grid.best_score_)\n",
    "\n",
    "SFS2 = SFS(iq_best, k_features=(1, 14), forward=True, floating=True,\n",
    "           verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)\n",
    "sfs2 = SFS2.fit(X_iq.as_matrix(), y_iq)\n",
    "print([X_iq.columns[i-1] for i in sfs2.k_features])\n",
    "print(sfs2.k_score_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.fit(X_iq_base, y_iq_base)\n",
    "iq_base_best = grid.best_estimator_\n",
    "pred_iq_base = iq_base_best.predict(test_iq_base)\n",
    "print(grid.best_score_)\n",
    "\n",
    "SFS3 = SFS(iq_base_best, k_features=(1, 14), forward=True, floating=True,\n",
    "           verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)\n",
    "sfs3 = SFS3.fit(X_iq_base.as_matrix(), y_iq_base)\n",
    "print([X_iq_base.columns[i-1] for i in sfs3.k_features])\n",
    "print(sfs3.k_score_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Мы получили интересный результат: кросс-валидация показала, что для обоих городов лучший набор признаков - это ['precipitation_amt_mm', 'year']."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Гипотеза о том, что несмещенные признаки хуже для Икитос не подтверждается."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим модель с новыми признаками и узнаем значение метрики для нее:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.fit(X_sj[['precipitation_amt_mm', 'year']], y_sj)\n",
    "pred_sj = grid.best_estimator_.predict(test_sj[['precipitation_amt_mm', 'year']])\n",
    "\n",
    "grid.fit(X_iq_base[['precipitation_amt_mm', 'year']], y_iq_base)\n",
    "pred_iq_base = grid.best_estimator_.predict(test_iq_base[['precipitation_amt_mm', 'year']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base1 = pd.read_csv('/media/sadworker/DataRed/data/dengai/submission.csv', index_col=[0, 1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base1.total_cases = np.concatenate([np.round(np.exp(pred_sj) - 1).astype(int), \n",
    "                                   np.round(np.exp(pred_iq_base) - 1).astype(int)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base1.to_csv('/media/sadworker/DataRed/data/dengai/subm1imp.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Попробуем использовать в качестве обучающей выборку из 4 наиболее коррелирующих с целевым признаков:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_best.fit(X_sj[['year', 'reanalysis_min_air_temp_k', \n",
    "              'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], y_sj)\n",
    "pred2_sj = sj_best.predict(test_sj[['year', 'reanalysis_min_air_temp_k', \n",
    "                                   'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iq_base_best.fit(X_iq_base[['year', 'reanalysis_min_air_temp_k', \n",
    "                    'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], y_iq_base)\n",
    "pred2_iq_base = iq_base_best.predict(test_iq_base[['year', 'reanalysis_min_air_temp_k', \n",
    "                                                   'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base2 = pd.read_csv('/media/sadworker/DataRed/data/dengai/submission.csv', index_col=[0, 1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base2.total_cases = np.concatenate([np.round(np.exp(pred2_sj) - 1).astype(int), \n",
    "                                   np.round(np.exp(pred2_iq_base) - 1).astype(int)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base2.to_csv('/media/sadworker/DataRed/data/dengai/subm1corr.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Повторим те же операции, но уже с моделью Лассо:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_params = {'alpha': np.logspace(-10, 4, 50)}\n",
    "clf_lasso = Lasso(random_state = 13)\n",
    "gridl = GridSearchCV(clf_lasso, \n",
    "                    set_params,\n",
    "                    cv=tss, \n",
    "                    scoring='neg_mean_absolute_error',\n",
    "                    n_jobs=-1\n",
    "                   )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridl.fit(X_sj, y_sj)\n",
    "predl_sj = gridl.best_estimator_.predict(test_sj)\n",
    "\n",
    "SFSl1 = SFS(gridl.best_estimator_, k_features=(1, 14), forward=True, floating=True,\n",
    "           verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)\n",
    "sfsl1 = SFSl1.fit(X_sj.as_matrix(), y_sj)\n",
    "print([X_sg.columns[i-1] for i in sfsl1.k_features])\n",
    "print(sfsl1.k_score_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridl.fit(X_iq, y_iq)\n",
    "predl_iq = gridl.best_estimator_.predict(test_iq)\n",
    "print(grid.best_score_)\n",
    "\n",
    "SFSl2 = SFS(gridl.best_estimator_, k_features=(1, 14), forward=True, floating=True,\n",
    "           verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)\n",
    "sfsl2 = SFSl2.fit(X_iq.as_matrix(), y_iq)\n",
    "print([X_iq.columns[i-1] for i in sfsl2.k_features])\n",
    "print(sfsl2.k_score_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridl.fit(X_iq_base, y_iq_base)\n",
    "predl_iq_base = gridl.best_estimator_.predict(test_iq_base)\n",
    "print(gridl.best_score_)\n",
    "\n",
    "SFSl3 = SFS(gridl.best_estimator_, k_features=(1, 14), forward=True, floating=True,\n",
    "           verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)\n",
    "sfsl3 = SFSl3.fit(X_iq_base.as_matrix(), y_iq_base)\n",
    "print([X_iq_base.columns[i-1] for i in sfsl3.k_features])\n",
    "print(sfsl3.k_score_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лассо уступает xgboost на кросс-валидации."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9. Прогноз для тестовой или отложенной выборки (2 балла)\n",
    "\n",
    "    (+) Указаны результаты на тестовой выборке или LB score. Результаты на тестовой выборке сравнимы с результатами на кросс-валидации. Если тестовая выборка создавалась автором проекта, то механизм создания должен быть непредвзят и объяснен (применен разумный механизм выборки, в простейшем случае – рандомизация);\n",
    "    \n",
    "    (+/-) Значения метрик на тестовой выборке не сильно отличаются от значений метрик на кросс-валидации и/или тестовая выборка создана предвзято, но есть разумное обоснование этому (пример: заказчик взял тестовую выборку из другого распределения);\n",
    "    \n",
    "    (-/+) Значения метрик на тестовой выборке сильно отличаются от значений метрик на кросс-валидации и/или тестовая выборка создана предвзято;\n",
    "    \n",
    "    (-) Прогноз для тестовой или отложенной выборки отсутствует."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Модель xgbregressor со всеми признаками показала наилучший результат: MAE=27.6514."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Модель с двумя признаками, отобранными с помощью mlxtend не дала улучшения: MAE=34.0433"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Модель с коррелирующими признаками дала: MAE=30.1058"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Модель Лассо отстает от xgb на 2-4 единицы по метрике по всем наборам признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на результаты по отложенной выборке:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sjx_valid_df = pd.DataFrame(sjx_train[vid_sj:id_sj], columns=features_to_scale)\n",
    "sjx_valid_df['y'] = sj_df.total_cases[vid_sj:id_sj].values\n",
    "sjx_valid_df['y_adj'] = sj_df.total_cases_adj[vid_sj:id_sj].values\n",
    "\n",
    "sjx_train_df = pd.DataFrame(sjx_train[:vid_sj], columns=features_to_scale)\n",
    "sjx_train_df['y'] = sj_df.total_cases[:vid_sj].values\n",
    "sjx_train_df['y_adj'] = sj_df.total_cases_adj[:vid_sj].values\n",
    "\n",
    "sjx_test_df = pd.DataFrame(sjx_test, columns=features_to_scale)\n",
    "\n",
    "iqx_valid_df = pd.DataFrame(iqx_train[vid_iq:id_iq], columns=features_to_scale)\n",
    "iqx_valid_df['y'] = iq_df.total_cases[vid_iq:id_iq].values\n",
    "iqx_valid_df['y_adj'] = iq_df.total_cases_adj[vid_iq:id_iq].values\n",
    "\n",
    "iqx_train_df = pd.DataFrame(iqx_train[:vid_iq], columns=features_to_scale)\n",
    "iqx_train_df['y'] = iq_df.total_cases[:vid_iq].values\n",
    "iqx_train_df['y_adj'] = iq_df.total_cases_adj[:vid_iq].values\n",
    "\n",
    "iqx_test_df = pd.DataFrame(iqx_test, columns=features_to_scale)\n",
    "\n",
    "iqx_base_valid_df = pd.DataFrame(iq_base_train[vid_iq:id_iq], columns=features_to_scale)\n",
    "iqx_base_valid_df['y'] = iq_base_df.total_cases[vid_iq:id_iq].values\n",
    "iqx_base_valid_df['y_adj'] = iq_base_df.total_cases_adj[vid_iq:id_iq].values\n",
    "\n",
    "iqx_base_train_df = pd.DataFrame(iq_base_train[:vid_iq], columns=features_to_scale)\n",
    "iqx_base_train_df['y'] = iq_base_df.total_cases[:vid_iq].values\n",
    "iqx_base_train_df['y_adj'] = iq_base_df.total_cases_adj[:vid_iq].values\n",
    "\n",
    "iqx_base_test_df = pd.DataFrame(iq_base_test, columns=features_to_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_best.fit(sjx_train_df.drop(['y', 'y_adj'], axis=1), sjx_train_df['y_adj'])\n",
    "valf_sj = sj_best.predict(sjx_valid_df.drop(['y', 'y_adj'], axis=1))\n",
    "print('МАЕ на отложенной выборке для всех признаков и Сан-Хуана равно:%0.3f' % mean_absolute_error(sjx_valid_df['y'],\n",
    "                                                                                                np.exp(valf_sj)-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iq_base_best.fit(iqx_base_train_df.drop(['y', 'y_adj'], axis=1), iqx_base_train_df['y_adj'])\n",
    "valf_iq_base = iq_base_best.predict(iqx_base_valid_df.drop(['y', 'y_adj'], axis=1))\n",
    "print('МАЕ на отложенной выборке для всех признаков и Икитос равно:%0.3f' % mean_absolute_error(iqx_base_valid_df['y'],\n",
    "                                                                                             np.exp(valf_iq_base)-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_best.fit(sjx_train_df[['precipitation_amt_mm', 'year']], sjx_train_df['y_adj'])\n",
    "val2_sj = sj_best.predict(sjx_valid_df[['precipitation_amt_mm', 'year']])\n",
    "print('МАЕ на отложенной выборке для 2 признаков и Сан-Хуана равно:%0.3f' % mean_absolute_error(sjx_valid_df['y'],\n",
    "                                                                                                np.exp(val2_sj)-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iq_base_best.fit(iqx_base_train_df[['precipitation_amt_mm', 'year']], iqx_base_train_df['y_adj'])\n",
    "val2_iq_base = iq_base_best.predict(iqx_base_valid_df[['precipitation_amt_mm', 'year']])\n",
    "print('МАЕ на отложенной выборке для 2 признаков и Икитос равно:%0.3f' % mean_absolute_error(iqx_base_valid_df['y'],\n",
    "                                                                               np.exp(val2_iq_base)-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sj_best.fit(sjx_train_df[['year', 'reanalysis_min_air_temp_k', \n",
    "                          'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], sjx_train_df['y_adj'])\n",
    "val4_sj = sj_best.predict(sjx_valid_df[['year', 'reanalysis_min_air_temp_k', \n",
    "                                        'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])\n",
    "print('МАЕ на отложенной выборке для 4 признаков и Сан-Хуана равно:%0.3f' % mean_absolute_error(sjx_valid_df['y'],\n",
    "                                                                                                np.exp(val4_sj)-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iq_base_best.fit(iqx_base_train_df[['year', 'reanalysis_min_air_temp_k', \n",
    "                                    'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], iqx_base_train_df['y_adj'])\n",
    "val4_iq_base = iq_base_best.predict(iqx_base_valid_df[['year', 'reanalysis_min_air_temp_k', \n",
    "                                                       'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])\n",
    "print('МАЕ на отложенной выборке для 4 признаков и Икитос равно:%0.3f' % mean_absolute_error(iqx_base_valid_df['y'],\n",
    "                                                                                             np.exp(val4_iq_base)-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Таким образом, результаты на отложенной выборке показывают плохую корреляцию с тестовыми данными."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В условиях малой выборки и шумных данных единственным способом повышать результат видится валидирование по Leaderboard."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10. Выводы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В рамках проекта исследованы предоставленные данные, выявлены отдельные тренды и использованы 2 модели с подбором гиперпараметров.\n",
    "\n",
    "На данный момент результат на Leaderboard очень посредственный - только топ-75%. В процессе кросс-валидации сделан вывод о ее сомнительной ценности для решения данной задачи."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В дальнейшем стоит:\n",
    "* продолжить исследование, связанное со сдвигом данных относительно целевой переменной, ведь для этого есть теоретическое обоснование\n",
    "* исследовать и попытаться внести в модель выявленную цикличность целевой переменной в течение года, а также межгодовые тренды. Возможно, просто применить статистические методы.\n",
    "* рассмотреть использование ансамбля моделей"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
